Zadania

Zadanie 1 (1p)

Porównać w języku Julia reprezentację bitową liczby 1/3 dla Float16, Float32, Float64 oraz liczby, która jest inicjalizowana jako Float16, a potem rzutowana na Float64.

a = Float16(1/3)
b = Float32(1/3)
c = Float64(1/3)
d = Float64(Float16(1/3))

println("Float16(1/3) = ", a, " ", bitstring(a))
println("Float32(1/3) = ", b, " ", bitstring(b))
println("Float64(1/3) = ", c, " ", bitstring(c))
println("Float64(Float16(1/3)) = ", d, " ", bitstring(d))

Float16(1/3) = 0.3333 0011010101010101
Float32(1/3) = 0.33333334 00111110101010101010101010101011
Float64(1/3) = 0.3333333333333333 0011111111010101010101010101010101010101010101010101010101010101
Float64(Float16(1/3)) = 0.333251953125 0011111111010101010101000000000000000000000000000000000000000000

Zadanie 2 (2p)

Zbadać, jak zmienia się odległość między kolejnymi liczbami zminnoprzecinkowymi reprezentowanymi w komputerze za pomocą języka Julia. Narysować wykres używając Plots zależności odległości od wartości liczby dla zakresu od 1.0 do 1000000.0.

using Plots

x = 1:100:1000000
dist32(x) = eps(Float32(x))
dist64(x) = eps(Float64(x))
y32 = dist32.(x)
y64 = dist64.(x)

scatter(x, [y32 y64], xlabel="x", ylabel="distance", title=["Float32(x)" "Float64(x)"], legend=false, layout=2, colour=[:red :blue])

Zadanie 3 (2p)

Jedną z bibliotek numerycznych, jaką dodatkowo będziemy używać na zajęciach jest GSL (język C). Opis jak używać . Korzystając ze wsparcia dla wyświetlania reprezentacji liczb zmiennoprzecinkowych zobaczyć jak zmienia się cecha i mantysa dla coraz mniejszych liczb. Zaobserwować, kiedy matysa przestaje być znormalizowana i dlaczego?

Kod załączyć jako komórka Markdown sformatowana jako C (link). Wynik także jako Markdown (kod albo fragment zrzutu ekranu).

#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>

int main (void)
{
    float x = 1e-32;

    while(x > 0) {

        x = x/(2.0);

        printf ("x = ");
        gsl_ieee_printf_float(&x);
        printf ("\n");
    }

    return 0;
}

x =  1.10011111011000100011111*2^-108
x =  1.10011111011000100011111*2^-109
x =  1.10011111011000100011111*2^-110
x =  1.10011111011000100011111*2^-111
x =  1.10011111011000100011111*2^-112
x =  1.10011111011000100011111*2^-113
x =  1.10011111011000100011111*2^-114
x =  1.10011111011000100011111*2^-115
x =  1.10011111011000100011111*2^-116
x =  1.10011111011000100011111*2^-117
x =  1.10011111011000100011111*2^-118
x =  1.10011111011000100011111*2^-119
x =  1.10011111011000100011111*2^-120
x =  1.10011111011000100011111*2^-121
x =  1.10011111011000100011111*2^-122
x =  1.10011111011000100011111*2^-123
x =  1.10011111011000100011111*2^-124
x =  1.10011111011000100011111*2^-125
x =  1.10011111011000100011111*2^-126
x =  0.11001111101100010010000*2^-126
x =  0.01100111110110001001000*2^-126
x =  0.00110011111011000100100*2^-126
x =  0.00011001111101100010010*2^-126
x =  0.00001100111110110001001*2^-126
x =  0.00000110011111011000100*2^-126
x =  0.00000011001111101100010*2^-126
x =  0.00000001100111110110001*2^-126
x =  0.00000000110011111011000*2^-126
x =  0.00000000011001111101100*2^-126
x =  0.00000000001100111110110*2^-126
x =  0.00000000000110011111011*2^-126
x =  0.00000000000011001111110*2^-126
x =  0.00000000000001100111111*2^-126
x =  0.00000000000000110100000*2^-126
x =  0.00000000000000011010000*2^-126
x =  0.00000000000000001101000*2^-126
x =  0.00000000000000000110100*2^-126
x =  0.00000000000000000011010*2^-126
x =  0.00000000000000000001101*2^-126
x =  0.00000000000000000000110*2^-126
x =  0.00000000000000000000011*2^-126
x =  0.00000000000000000000010*2^-126
x =  0.00000000000000000000001*2^-126
x =  0

// mantysa przestaje być znormalizowana po przekroczeniu zakresu cechy

Zadanie 4 (4p)

Na przykładzie wybranego algorytmu niestabilnego numerycznie:

1. Pokazać, że działa źle.
2. Pokazać które konkretnie działania powodują zwiększenie błędu (np. dzielenie przez małą liczbę, cancellation).
3. Zademonstować wersję stabilną.

Nie wolno pokazywać przykładów z wykładu (lub bardzo podobnych)!

Wszystkie punkty przedstawić w postaci notatnika Julii.

# Funkcja:
f(x) = cos(x) ^ 2

# Algorytm niestabilny numerycznie:
# Funkcja obliczająca pochodną funkcji

function f_derivative(f, x, e)
    return (f(x) - f(x-e)) / e
end

# 1:
println("Oczekiwany wynik: -2 * cos(1) * sin(1) ~= ")
println(-2 * cos(1) * sin(1))
println()

println("Wyniki algorytmu niestabilnego:")
println(f_derivative(f, 1, 0.001))
println(f_derivative(f, 1, 0.0000001))
println(f_derivative(f, 1, 0.0000000000001))
println(f_derivative(f, 1, 0.00000000000001))
println(f_derivative(f, 1, 0.000000000000001))
println(f_derivative(f, 1, 0.0000000000000001))
println()

# 2. Dobranie niewłaściwej wartości epsilon (e) powoduje zwiększanie błędu.
# Zbyt duże e, rzędu ~0.001 powoduje odbiegający od rzeczywistości wynik algorytmu, a zbyt małe e mniejsze od systemowego epsilon, może spowodować "wyzerowanie" wyniku.

# 3. Wersja stabilna:

function f_derivative2(f, x)
    e = sqrt(eps(Float64))
    return (f(x+e) - f(x-e)) / (2e)
end

println("Wynik algorytmu stabilnego:")
println(f_derivative2(f, 1))

Oczekiwany wynik: -2 * cos(1) * sin(1) ~= 
-0.9092974268256818

Wyniki algorytmu niestabilnego:
-0.9097129673252824
-0.909297467877046
-0.9092726571680032
-0.8992806499463768
-0.8326672684688673
0.0

Wynik algorytmu stabilnego:
-0.9092974290251732